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Abstract 

Quantum Optimal Control Theory (QOCT) provides the necessary tools to theoretically design driving 
fields capable of controlling a quantum system towards a given state or along a prescribed path in Hilbert 
space. This theory must be complemented with a suitable model for describing the dynamics of the quan- 
tum system. Here, we are concerned with many electron systems (atoms, molecules, quantum dots, etc) 
irradiated with laser pulses. The full solution of the many electron Schrodinger equation is not feasible in 
general, and therefore, if we aim to an ab initio description, a suitable choice is time-dependent density- 
functional theory (TDDFT). In this work, we establish the equations that combine TDDFT with QOCT, and 
demonstrate their numerical feasibility with examples. 

PACS numbers: 42.50.Ct, 32.80.Qk, 02.60.Pn 
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The quest for systems able to perform quantum computing yj], the synthesis of design- 
molecules by laser-induced chemical reactions 121], or the control of electron currents in molecular 
switches using light [3] may benefit from the recent advances in the field of design and synthesis of 
laser pulses specially tailored to perform a specific task [4]. The laser pulse creation and shaping 
techniques have improved impressively over the last decades, and thus the area of experimental 
optimal control has become a well established field. 

Suchpulses can also be theoretically derived with the help of quantum optimal control the- 
ory |5kZ|] (QOCT). This theory is rather general in scope, and its basic formulation makes no 
assumption on the nature and modelling of the quantum system on which the pulse is applied. In 
practice, the solution of the QOCT equations requires multiple propagations, both forwards and 
backwards, for the system under study. Since these propagations are in general unfeasible for 
many-particle systems, few-level simplifications and models are typically postulated when han- 
dling the QOCT equations. Unfortunately, these simplifications are not always accurate enough: 
strong pulses naturally involve many levels, and normally perturbative treatments are not useful. 
Non linear laser-matter interaction must sometimes be described ab initio. 

In this work, we are concerned with many-electrons systems (atoms, molecules, quantum 
dots...) irradiated with femtosecond (or attosecond) pulses, with intensities typically ranging 
from 10 11 to 10 15 Wcm" 2 (a non-linear regime that nevertheless allows for a non-relativistic treat- 
ment). This may lead to a number of interesting phenomena, e.g. above-threshold or tunnel 
ionization, bond hardening or softening, high harmonic generation, photo-isomerization, photo- 
fragmentation, Coulomb explosion, etc [8]. In order to describe this type of processes from first 
principles, time-dependent density-functional theory [9] (TDDFT) has emerged as a viable alter- 
native to more computationally expensive approaches based on the wave function. 

In TDDFT, the system of interacting electrons is substituted by a proxy system of non- 
interacting electrons - the "Kohn-Sham" (KS) system, which is computationally much less de- 
manding. The theory guarantees the identity of the electronic densities of the two systems, and 
the existence of a density functional for each possible observable, thus allowing (in principle) the 
computation of any property without ever dealing with the many-body wave function. The theory 
is however hindered by the lack of knowledge of the precise external potential seen by the auxiliary 
non-interacting system (the so-called "exchange and correlation" potential, which is a functional 
of the density itself, has to be approximated), and, in many cases, the precise form of the density 
functional that provides the required observable. Fortunately, a number of valid approximations 



for these density functionals have been developed over the years, which have made of TDDFT a 
computationally efficient possibility to describe many processes. 

We are thus led to the necessity of inscribing TDDFT into the general QOCT framework. We 
will lay down and discuss the equations that result when TDDFT is used to model the system. 
Then, in order to demonstrate the computational feasibility, we present one sample calculation: a 
2D two-electron system being optimally driven between two potential wells. 

In the spirit of TDDFT, we substitute the problem of formulating QOCT in terms of the real 
interacting system, by the formulation of the optimization problem for the non-interacting system 
of electrons. The equations of motion for the single-particle orbitals of this system, also known as 
time-dependent KS (TDKS) equations, are: 

i^(ro,0 = (ro\H KS [n x(i} (t),u,t}\<Pi(t)) = 
- 1 V 2 (p; (rc,t) + (re | % | % (t)) + v H [n(t)} (?) q>,- (rc,t)+ 

(ro\V xc [n za (t)]\q>i(t)) + (ro\V ext [u]\<pi(t)) , (1) 

N 

n x(a (r,t) = £9*(rT,0cp;(**V) , (2) 
(=l 

n(f,t) = J^naa(f,t) , (3) 
a 

for i = 1 , . . . , N orbitals. The greek indexes o, x, (0 . . . run over the two spin configurations, up and 
down. The densities are, by construction, equal to that of the real, interacting system of electrons. 
% represents the internal, time independent fields - usually a nuclear Coulomb potential V n (r), 
and may include as well a spin-orbit coupling term of the form o • W n x p (where a is the vector 
of Pauli matrices). The term v#[n(f)](r) = JdV "i^j is the Hartree potential, and Vxcfntco] is the 
exchange and correlation potential operator, whose action is given by: 

(?a|y xc [M0]|(p ( -(0> = £vS[n T(n (0](?)cpK^,0 • (4) 

P 

Note that, for the sake of generality, we allow a spin-resolved exchange and correlation potential, 
that depends on the four spin components n T § (in many situations more restricted dependences are 
assumed). However, we do assume here an adiabatic approximation, i.e. v xc at each time t is a 
functional of the densities at that time, n X (o(t). This restriction is non-essential for the derivations 
that follow, but the use of non-adiabatic functionals is very scarce, and adiabatic approximations 
will result in simpler equations. 
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The last potential term, V ext , is the external time-dependent potential, which is determined by a 
"control" u. In a typical case, this external potential is the electric pulse created by a laser source 
in the dipole approximation, and u is the real time-dependent function that determines its temporal 
shape (in this case, (Fo|y e xt[w] 19/(0 ) = u {t)r- Ptyiiroj), where p is the polarization vector of the 
pulse). We write it however in general operator form, since it can be a 2x2 matrix which may 
include both a time-dependent electric field as well as a Zeeman-coupled magnetic field [Q]. The 
mathematical nature of u may also be diverse: it may not be a time-dependent function, but a set 
of N parameters that determine the precise form of the electric field. 

If we group the N single particle states into a vector cp(f), we can rewrite the TDKS equations 
in a matrix form: 

i9(0 = g[n™(0,w,*]<f>(0) (5) 

where H_[n Z(£l (t) , u, t] = #ks [ n tco(0 5 u i ^Ly an ^ L N is the A^-dimensional unit matrix. With this nota- 
tion we stress the fact that we have only one dynamical system - and not N independent ones, since 
all cp ; are coupled. This coupling, however, comes solely through the density, since the Hamilton 
matrix is diagonal. 

The specification of the value of the control u, together with the initial conditions, determines 
the solution orbitals: u — >■ <p[w]. Our task is now the following: we wish to find an external field - in 
the language of OCT, a control u - that induces some given behaviour of the system, which can be 
mathematically formulated by stating that the induced dynamics maximizes some target functional 
F . Since we are using TDDFT, this functional will be defined in terms of the KS orbitals, and will 
possibly depend explicitly on the control u: 

F = F[<p,u]. (6) 

Since the KS orbitals depend on u as well, the goal of QOCT can be formulated as finding the 
maximum of the function: 

G[u] = F[(p[u],u]. (7) 

In the most general case, the functional F depends on cp at all times during the process (we have a 
"time-dependent target"). In many cases, however, the goal is the achievement of some target at a 
given time T that determines the end of the propagation interval (we then have a "static target"). 
In both cases, the determination of the value of the function G is obtained by performing the 
propagation of the system with the field determined by the control u. 
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There are many optimization algorithms capable of maximizing functions utilizing solely the 
knowledge of the function values ("gradient-free algorithms"). We have recently employed one of 
them in this context [11]. However, QOCT provides the solution to the problem of computing the 
gradient of G - or, properly speaking, the functional derivative if u is a function. The non-linear 
dependence of the Hamiltonian with the density slightly complicates the derivation, but we sketch 
the key steps: First, we must note that searching for a maximum of G is equivalent to a constrained 
search for F - constrained by the fact that the cp orbitals must fulfill the TDKS equations. In order 
to do so, we introduce a new set of orbitals % that act as Lagrange multipliers, and define a new 
functional 7 by adding a Lagrangian term L to F: 



J[q>,%,u] =F[q>,«]+L[<p,x,M] 



T d 

d* (Xj(t) |i— + H KS [nuu(t) , u, t]\<$>j(t) ) 



(8) 
(9) 



N 

= -2£Re 

7=1 

Setting the functional derivatives of 7 with respect to the % orbitals to zero, we retrieve the TDKS 
equations. In an analogous manner, we obtain a set of solution %[u] orbitals by taking functional 
derivatives with respect to op: 
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X(T) = 0. 

The non-diagonal operator matrix |c [<p[w] (t)] is defined by the operators: 

(ro|ly[<p[ M ](?)]|\|/) = -2i^cp i -[ M ](rK,?): 



(10) 
(11) 
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If we now note that G[u] = J[(p[u],(p[u},u], we arrive to: 



(12) 



(13) 



V u G[u] =V B F[<p,u] 
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Several aspects of these equations deserve further discussion: 

(1) Eqs. [10] and [TTJare a set of first-order differential equations, whose solution must be ob- 
tained by backwards propagation, since the boundary condition, Eq.QjIis given at the end of the 
propagating interval, T. Note that this propagation depends on the time-dependent KS orbitals 
(p[w]. Therefore, the numerical procedure to follow consists of a forward propagation to obtain 
cp[w], followed by a backwards propagation to obtain %[u]. 

(2) These backwards equations are non-homogeneous, due to the presence of the last term in 
Eq. [RE the functional derivative of F with respect to cp - but see point (4) below. 

(3) Often, the control target functional F is split like: 

F[q>,K]=Ji[<p]+J 2 [«]. (15) 

J\ codifies the actual purpose of the optimization, whereas Ji imposes a penalty on the control 
function, in order to avoid, for example, the solution field to have unreasonable amplitudes. In the 
following, we will assume this division. 

(4) The previous equations (110H1 II) refer to a general "time-dependet target" case, as mentioned 
above. In many cases of interest, the target functional F takes a "static" form, which can be 
expressed as: 

J 1 [<?,u}=0[q(T),u], (16) 

for some functional O whose argument is not the full evolution of the KS system, but only its value 
at the end of the propagation. In this case, the non-inhomogeneity in Eq. [10] vanishes, and instead 
we obtain a different final-value condition: 

^'^ = -^r- <17) 

(5) The previous Eq. [14] assumes that u is a set of N parameters, u G R that determines the 
control function. If u is directly the control function, the gradient has to be substituted by a 
functional derivative, and the result is: 



5G 5F[<p,n] 



8u(t) 8u(t) 



+ 2Im 

<p=<pM 



£(x>](0|£|cp>](0> 

.7=1 



(18) 



We have assumed here that the external potential v ext is determined by the function u by a linear 
relationship: 

V ext [u](t) = u{t)D. (19) 



This is the most usual case would be the dipole operator, and u(t) the amplitude of an electric 
field), but of course it would be trivial to generalize this to other possibilities. 

The previous scheme permits therefore to control the KS system. However, the goal is to 
control the real system. In principle, the target is given by some functional J\ [*P] that depends 
on the real many-electron wave function of the interacting system. This object is not provided by 
TDDFT, that only provides the density n. Therefore, the ideal situation would be that in which J\ 
depends on *F only through the density n, J\ = J\ [n]. In this manner, optimizing for the KS system 
is strictly equivalent to optimizing for the real one. For example, this holds if J\ is given by the 
expectation value of some one-body local operator A: 

J\[y] = (y(T)\A\V(T)) = J d 3 rn(r,T)a(r), (20) 

where A = m tn i s case, Eq.[T7]is simply: 

X i [u]{fa,T)=a{f)a? i [u]{raJ). (21) 

This will be the kind of target that we will be using in the example below. Note that TDDFT 
ensures that all observables are functionals of the density, and therefore in principle one could 
always write such functionals (we will provide another example in a forthcoming publication, in 
which the target is the high harmonic generation spectrum, which is also an explicit functional of 
the density). 

Unfortunately, in some cases we still need to find the explicit density functional dependence 
in many cases of interest. For example, a very common control goal is the transition from an 
initial state to a target state. In other words, the control operator A is the projection operator onto 
the target state A = ^target) (^target |- We have no exact manner to substitute, in this case, the J\ 
functional by a functional J\ defined in terms of the density, or in terms of the KS determinant. It 
can be approximated, however, by an expression in the form: 

JM = mT)\^ Cl W)\\ (22) 
/ 

where cp(T) is the TDKS determinant at time T, and we compute its overlap with a linear com- 
bination of Slater determinants cp 7 , weighted with some coefficients cj. These Slater determinants 
would be composed of occupied and unoccupied ground state KS orbitals, cp 7 = det[cp 7 , . . . ,<$L]. 
In this case, Eq.[l7]takes the form: 

Xi[u](ro,T) =£^(ro)((p(r)|( P / )(cp 7 |9(n), (23) 
u 
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FIG. 1. (a) External potential defining a model for double quantum dot. (b) Density of the initial, ground 
state (blue) and final, propagated density (red), (c) Optimized electric field for the charge-transfer process 
described in the text, (d) Convergence history of the conjugate gradient algorithm. All magnitudes are given 
in effective atomic units. 

\ u (r<5) = c^TrUM'r'A* (Fo)} , (24) 

where M' mn = (<p m |q>£) and A}(ro) m „ = 8 mi (^ n (fo). 

We have implemented the described TDDFT+QOCT formalims in the octopus code |l2|1 . In 
the following, we describe a simple example: the charge transfer between two neighboring poten- 
tial wells, considered as models for 2D quantum dots, such as the ones created in semiconductor 
heterostructures. We consider a two-electron system, trapped in an asymmetric double quantum 
dot well modelled by a potential function given by (in the following, we consider effective atomic 
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units): 

v ^,y) = ^x 4 - 1 -x 2 + ^+ 1 -y 2 . (25) 

The potential landscape is depicted in Fig.QJa). We then solve the ground state KS equations for 
this system, by making use of the local density approximation to the exchange and correlation 
parameterized by Attacalite et al ft 1 311 . The ground state density will be localized in the left well 
(see Fig.QJb): n gs (r) = 2 1 cp§ s (F) | 2 , where (p[j s is the ground state KS orbital. 

We apply an electric field, polarized along the x direction. Its amplitude is parameterized by its 
Fourier coefficients: 



M/2 



,2 fin \ , 2 . (2% 
a„\ —cos — nt + b n \ — sin — nt 



(26) 



The {a n ,b n } coefficients are therefore the control u (although the constraint Y^j2\ a n = is en- 
forced in order to ensure £(0) = e(T) = 0). Since our goal is to transfer as much charge as possible 
from the left to the right well, we formulate a target in the form: 

M/2 

rn(r,T) — 
'*>° n=\ 

In words, we wish to arrive to a state in which all the density is localized in the x > region. The 
last term of Eq. ITTI corresponds to the penalty: 

,T M/2 



M/2 

F[y\a,b]= / d 3 rn(r,r)-a£ (a 2 n + b 2 n ) . (27) 

- Jx>0 



h\a,b\ = -a f dte 2 (t) = -a £ (a% + b%) 

J° n=l 



(28) 



and it is introduced in order to prevent the solution field from having too much intensity. 

The solution field is shown in Fig. Etc). We have employed a standard conjugate gradients 



(CG) algorithm to perform the optimization. After around 60 CG iterations [14], the control field 
is converged and we achieve a value of 1.92 for 7i - the maximum is 2 (see convergence plot in 

Kg. ID- 

In conclusion, we have shown how TDDFT can be combined with QOCT, and we have demon- 
strated how the resulting equations are numerically tractable. This provides a scheme to perform 
QOCT calculations from first principles, in order to obtain tailored function- specific laser pulses 
capable of controlling the electronic state of atoms, molecules, or quantum dots. Most of the pre- 
vious applications of QOCT were targeted to control, with femto-second pulses, the motion of 
the nuclear wave packet on one or few potential energy surfaces, which typically happens on a 
time scale of hundreds of femtoseconds or picoseconds. The approach developed in this Letter, 



on the other hand, is particularly suited to control the motion of the electronic degrees of freedom 
which is governed by the sub-femto- second time scale. The possibilities that are open thanks to 
this technique are numerous: shaping of the high harmonic generation spectrum (i.e. quenching 
or increasing given harmonic orders), selective excitation of electronic excited states that are oth- 
erwise difficult to reach with conventional pulses, control of the electronic current in molecular 
junctions, etc. Work along these lines is in progress. 

This work was partially supported by the Deutsche Forschungsgemeinschaft within the SFB 
658. 
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